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ABSTRACT 

A  finite  difference  approximation  to  a  non-linear 
set  of  parabolic  differential  equations  arising  in  shallow 
water  theory  is  given.   These  difference  equations  were 
used  to  determine  the  shape  and  rate  of  propagation  of  a 
hump  of  fluid  down  a  channel  of  constant  depth.   The  hump 
of  fluid  was  found  to  spread  instead  of  steepen,  as  is  the 
case  in  the  usual  shallow  water  theory. 
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1,   Introduction.   As  is  well  known  the  shallow  water  equations 
for  the  flow  of  fluid  down  a  channel  Is  hyperbolic  [5]  and  the 
theory  suffers  from  the  defect  that  any  hump  of  fluid  traveling 
down  a  channel  will  eventually  fall  over  because  the  top  of  the 
hump  moves  relatively  faster  than  the  trougho   Some  years  ago, 
Boussinesq  [l]  suggested  an  improvement  of  the  shallow  water 
theory  which  is  obtained  by  taking  greater  account  of  the  verti- 
cal component  of  the  velocity„   In  this  improvement,  Boussinesq 
arrived  at  the  following  system  of  non-hyperbolic  equations  to 
describe  the  flow  of  fluid  in  a  channels 

(1.1)  d     (|;^h)  ^  dj^^llhl  ^  0        (conservation) 

^^•2)  §T  -^  ^§7  -^  g  ^'(3x   •"  3  ^    2.   "  ^   (momentum). 

We  take,  as  usual,  ^  to  be  the  height  of  the  fluid  measured  from 
the  undisturbed  level   h  ,   u   to  be  the  horizontal  velocity  of 
the  fluid,   g   to  be  the  gravitational  constant,   x   to  be  the 
distance  downstream  and   t   to  be  the  time. 

Fig.  1 


This  study  was  suggested  by  Professor  J,  Jo  Stoker  who  was 
interested  in  determining  whether  the  equations  (1.1)  yield 
a  realistic  description  of  the  breaking  of  a  wave.   The  authors 
are  indebted  to  Professor  J,  J,  Stoker  for  his  advice  and 
encouragement  in  carrying  out  this  investigationo   Our 
calculations  indicate  that  no  waves  break  under  (lol)  in  contrast 
to  the  shallow  water  theory  (see  [5])  i^i  which  all  waves  break. 
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The  conservation  equation  (1.1)  is  identical  with  the 
conservation  equation  of  the  shallow  water  theory  [Ij  but 
the  momentum  equation  has  an  added  term  involving  a  third 
derivative  of   (i+h)  ,   Recently,  equations  with  a  third 
derivative  term  have  appeared  in  problems  of  magneto- 
hydrodynamics  [3].   Because  of  the  interest  in  equations  of 
this  "type"  and  the  lack  of  analytic  solutions  of  these  non- 
linear equations  it  is  desirable  to  develop  a  numerical 
technique  for  the  solution  of  such  equations » 

In  this  paper  we  discuss  an  implicit  numerical  scheme  for 
the  solution  of  equations  of  the  type  (l„l-2)<,   After  trans- 
forming the  equations  into  dimensionless  form,  we  change  the 
variables  and  recombine  the  equations  in  a  manner  that  would 
lead  to  the  characteristic  form  if  the  third  derivative  term 
were  not  present  in  the  momentum  equationo   The  initial  and 
boundary  conditions  assumed  to  be  appropriate  for  this  system 
are  those  required  by  the  shallow  water  theory. 

Experiments  which  were  performed  on  an  IBM  70L|.  digital 
computer  for  an  initial  boundary=value  problem  indicated  that 
numerical  errors  grew  from  the  left  boundary,  but  that  the 
scheme  is  stable  away  from  the  boundary.   In  order  to  over- 
come this  instability  we  modified  the  difference  equation  at 
the  boundary  using  the  conservation  equation  to  correct  the 
height  in  an  implicit  fashion,  i.e.  at  the  advanced  time  step. 
The  instability  from  the  left  boundary  then  occurred  one 
point  to  the  right  i.e„  by  using  the  conservation  ^f  mema 


equation  at  the  left  boundary  we  had  moved  the  cause  of  the 
perturbation  to  the  right.   We  then  used  the  implicitly 
differenced  conservation  of  mass  equation  at  several  points 
adjacent  to  the  left  boundary.   The  start  of  the  instability 
moved  correspondingly  to  the  right.   We  finally  eliminated 
the  instability  by  using  this  conservation  equation  for  all 
of  the  points o   That  is,  we  used  an  explicit  scheme  to  get  a 
first  estimate  of  the  unknown,  but  then  we  corrected  this 
estimate  with  an  implicit  scheme » 

With  this  method  some  numerical  solutions  of  the 
propagation  of  a  hump  of  fluid  into  a  channel  of  constant  depth 
were  obtained  on  the  IBM  TOi^..   These  solutions  indicate  that 
the  hump  of  fluid  does  not  steepen  but  spreads  in  a  manner 
found  in  an  analytic  approximation  by  Gardner  and  Morikawa  [3] 
for  a  similar  problem, 

2.   The  Equations  in  "characteristic"  form. 

If  the  third  derivative  term  were  missing  from  the 
momentum  equation  (lo2),  we  would  have  the  usual  hyperbolic 
system  of  shallow  water  theory.   Since  the  third  derivative 
term  arises  from  an  attempt  to  make  a  small  correction  to  the 
shallow  water  theory,  and  since  the  characteristic  form  is  very 
useful  in  the  discussion  of  the  shallow  water  theory,  we 
will  cast  the  equations  (1.1-2)  into  a  form  which  would  be 
the  characteristic  form  if  the  third  derivative  term  was  not 
present. 

First  let  us  transform  to  dlmensionless  variables  by  the 
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following  transformations 

(2.1)  K  ^P^'>   -^  ^ 

(2.2) u  — >  u 

\/gh 

(2.3)  vT^ 


With  these   transformations   equations    (1.1=2)    become 

(2.5)  ft*^     -      ° 

(2.6)  Jt        ^  fe       Jx       ~ ^     ~  • 

<^t  cjx       dx      ^^^^ 

By  using  the  conservation  equation  (2.5)  we  can  change  the 

third  derivative  term  to   -  V  C^^v   ,  a  third  derivative  term 

C)xcjt 

involving  only  one  time  differentiation. 

Wo  now  let 

(2.7)  c^  =  h 

and  obtain  from  {2,S-^)    a  form  similar  to  the  characteristic 
form  [5] 

(2.8)  J^.,u.c)^^j  (u.20).d2|^  =   0 

Even  though  the  third  derivative  term  may  physically  give 
a  small  correction  to  the  usual  shallow  water  equations,  a 
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possible  stable  difference  approximation  will  have  to  treat 
the  third  derivative  term  as  not  being  small,  and  hence  we 
will  have  to  use  the  third  derivative  term  in  solving  for 
the  unknown  quantities  at  the  advanced  time  step. 

Since  both  variables  appear  naturally  combined  in  the 
third  derivative  term,  we  define 
(2,10)  q  =  u  c^ 

called  the  discharge  (q  is  a  measure  of  the  volume  of  water 
passing  a  given  point  per-unit  of  time)  .   We  then  have 

=  0 


tit 

=  0 

If  the  third  derivative  term  were  absent,  we  would  have 
to  prescribe  the  appropriate  initial  conditions  for  a  solution 
in  the  form 

(2.13)    (loCe)       q(x,0)  =  given  ,  x(x,0)  =  given  ; 
and  if  we  had  a  semi-infinite  channel,  0  <  x  <  oo  we  would 
need  only  specify  one  of  the  two  unknowns  on  the  boundary 
X  =  0  (under  normal  circumstances  in  which  u+c  >  0  ,  u-c  <  0 
at  X  =0  cf.  [5])'   We  find  it  convenient  to  prescribe  the 
discharge  at  x  =  0 
(2.11^)     (B.C)        q(0,t)  =  given   . 

Since  we  ha/e  a  third  derivative  term  we  might  need 
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additional  conditions.   In  the  form  given  by  (2.11-12)  the 
spatial  derivative  is  raised  by  one  order  over  that  of  the  usual 
shallow  water  theory,  while  the  time  derivative  is  still  of  first 
order.   This  observation  suggests  that  only  an  extra  boundary 
condition  might  be  required.   It  seems  reasonable  to  take  the 
extra  boundary  condition  at  x  =  oo  .   We  will  be  interested  in 
creating  a  wave  motion  at  x  =  0  and  in  observing  the  motion 
of  the  wave  towards   x  =  oo  .   We  will  not  wish  any  disturbances 
to  arise  from  x  =  oo  and  hence  we  will  assume  that  for  fixed 
t  >  0  , 

(2.15)   (B.C.)   lim  q(x,t)  =  q(oo,  0)  ,   and   lim  q^(x,t)  =  0  , 
x-*»oo  x->oo 

3,   The  Finite  Difference  Scheme 

We  will  now  consider  a  finite  difference  approximation  to 
the  "characteristic  type"  equations  (2,11-12)  which  needs  only 
the  initial  and  boundary  conditions  given  by  (2.13-15) •   Two 
considerations  motivate  our  choice  of  difference  approximation: 

(I)  if  the  third  derivative  term  were  not  present,  the 
difference  approximation  should  be  a  known  stable 
approximation  to  the  shallow  water  theory  equations, 

(II)  the  third  derivative  term  should  be  used  appropriately  to 
advance  the   q  , 

We  consider  a  mesh  such  that  the  spatial  and  time  mesh  points 
are  designated  as   x,  ,  j  =  0,1,2,  ,..,  t^  ,  n  =  0,1,2,..,,   with 
Xq  =  0  as  the  left  boundary  and   t^  =  0  as  the  initial  time. 
The  unknowns   q(x, t)   and   c(x, t)   at  the  mesh  points  will  be 
written  with  capitals  and  subscripts: 
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(3.1)  ^(^j-V=V 

(3.2)  c(Xj,t^)  =  Cj^   . 

In  order  to  have  the  difference  approximations  of  (2.11-12) 

reduce  to  the  shallow  water  difference  approximations  when  the 

third  derivative  tenn  is  absent  we  will  take  backward  space 

differences  for  the  £^   (-^  +  2c)   term  in  (2.11)  and  forward 

space  differences  for  the  -^  (-^  -   2c)   term  in  (2.12).   This 

^        c 
difference  approximation  for  the  shallow  water  equations  is  known 

to  be  stable  (since   u+c  >  0  ,  u-c  <  0)  cf.  [2]).   To  treat  the 

third  derivative  term,  we  difference  it  forward  in  time  and  in 

a  centered  fashion  in  the  spatial  direction. 

Following  out  this  difference  approximation,  we  obtain  from 

(2.11-12) 

J  n+l  jn 

+   (__Ia  +  C  .  )  .-i-  (-4^  +  2  c ,     -  -4^^  -  2  C  .  ,  ) 
'q2  jn'Ax  q2       jn   ^2  j-ln^ 

jn  jn  j-ln 

^     (Q..n  „^n  "  2  Q,„^,  +  Q, 


(A^)  At 


^At        J+1  ^"*"1      jn+1   ^j-1  n+1 


-  S>ln  ^  2  «j„  -  «j.i„)  =  C  . 


Jn  jn 

Q  Q  Q 

q2  jn'  4x  V2         j+ln   ^2       jn 
Jn  j+ln  jn 


(^x)'l^t 


2    (^+1  n+1  -  2  Qjn+l  ^  S-1  n+1 


-  Q.  ,   +  2  Q.   -  Q,  -,  )  =  0 
"^j  +  ln   '^   "^jn   ^j-ln' 
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where  Ax  =  x  ^^  -  x^.   and  ^t   =  t^_^^  "  ^n  *   ^®  ^^^■'-  ^°-^^  ^^ 
and  At      constant  during  a  given  calculation.   The  equations 
(3.3-U-)  hold  for  all   j  >  1   and  all   n  >  0  .   For  n  =  0  we 
know  Q,.^  and   C  ,_   from  the  initial  conditions  (2.13).   For 
J  =  0  we  know   Q^.    from  the  boundary  condition  (2.1L|-),  and 
for   j   (  =  Jniax^   taking  on  its  maximum  value   j  =  j^^^  we 

know  Q.      from  the  boundary  condition  (2.15).   That  is,  we 

•'max 
approximate  the  infinite  channel  by  a  suitably  long  channel 

(the  length  of  the  channel  increases  as  the  calculation  proceeds). 

Upon  subtracting  (3.1|-)  from  (3.3),  we  easily  solve  for 


C  .   n   in  the  form: 
j  n+1 

^^'^"^      ^j  n+1  =  ^jn  "■  JJ^ 


Q 


i+ln 


jn   |_  j  +  ln 


jn     jn 


-  2(°J.l,n  -  ^J-l.n' 


+  c 


jn 


S-Hn   S-ln 


-7^^ 


^j  +  ln   ^j-m 


+  2  C  ,_^,  -2C  .   +  C  . 
j+ln   Jn    J 


..) 


=  0 


By  adding  (3.3-14-),  we  find  after  some  rearranging. 


(3.6) 


+  (. 


J  n+1 


r)  Q. 


Q. 


C:  _,.   (4x)2^  j  '^-^i   (4x)2  "J-i  "+i 


Sn   At   ^ 
Jn 


(^l±ln  .  ^jzlii)  +  2(C 


jn  ^    j+ln    j-ln 
+  C 


j+ln   jn   j-ln' 


jn 


(S-Hn     ^jn  ,  ^j-ln,   , 
c2        ^   c2    ^  '^^^j+ln^j-ln^ 
^^j  +  ln     ^jn   ^j-ln 


^  ^  (^j+ln  -  2  Qj,  .  Qj,,^)  , 


the  equations  for  Q,  ^^^    .   In  equations  (3.6)  the  unknowns 
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Q.      (for   j  =  1,  2s,  . . .,  ( J*  „„-l) )   appear  in  triplets  at  the 
J  n  I  J.  rna.A. 

advanced  time   t   ,  .   The  solution  of  this  tridlagonal  system 
can  be  easily  given  by  well-known  and  well-conditioned  recursion 
formulae. 

We  need  one  more  equation  to  complete  the  difference  scheme: 
an  equation  to  determine   C-   ,   for  we  have  assumed  that  only 
Qq   ,   is  given  on  the  left  boundary  and  the  difference  equation 
(3,5)  can  be  solved  for  the  unknowns   C.    .   for   i  =  1,2,  ...oi 
In  order  to  determine   C^   _^,   we  use  the  conservation  equation 
(2.5)  with  the  notation  (2.7)  and  take  a  one  sided  space  difference, 
to  find 

r  2     t  )  1/2 

^^-"^^  ^0  n+1  =   ^On  \^   ^^In  "  ^OnM 

We  recapitulate  our  algorithm  for  solving  the  difference 
equations  (3o5-7)  as  follows s 

(i)   We  use  (3„5)  to  obtain   C  ^^^  (j  >  1,  n  >  0)   (which 

can  be  done  since  all  the  quantities  on  the  right  side  of  (3.5) 

are  presumably  known  at  the  previous  time   t  )  .   Since  the 

right  boundary  is  at  infinity,  we  solve  for  C .   ,   for   j 

increasing  until   C.  .^  =  C.    to  within  some  preset  tolerance 
^         jn+1    jn  ^ 

and  define   i      to  be  this  last  value  of   1  » 

(ii)   We  use  (3o7)  to  obtain  C_   ,  . 

(lii)  We  solve  the  coupled  system  of  equations  (3o6)  for 

Q.   -,  6   The  boundary  conditions  on  the  left  and  right  are  just 

enough  for  the  solution  of  the  coupled  equations  (3.6),   For 

the  semi-infinite  channel  the  right  boundary  is  considered  to  be 

that  point  for  which   C.     =  C.    to  within  the  accuracy  of 

jn+1    jn 
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the  calculation. 

Numerical  calculations  for  large  discharges  at  the  left 
boundary  indicated  an  instability  near  the  left  boundary  and  no 
instability  away  from  the  left  boundary.   Since  the  instability 
showed  up  in  C  .   with  negligible  effect  on  Q,.   ,  we  were  led 
to  a  further  implicit  type  correction  for  C    „   The  implicit 
type  correction  calculation  for  C    that  we  employ  uses  the 
values  of  the  discharge  at  both  the  advanced  and  present  time 
step.   We  accomplish  the  desired  implicitness  by  differencing 
the  conservation  equation  (2,5)  as  follows; 

'3-8)   CjnH.1  =U%  -   ^  (9j.l  n.l  -  '^i-l   n.l  *  «J+ln  "  «J-ln' 
For  the  boundary  we  take  a  one  sided  space  difference: 

^•9)    ^On^l  =  f^On  -  ^   ^^m^l  -  V+1  "  ^in  '  %n^J   '  ^  ' 

Our  procedure  is  as  follows:   (i)  we  solve  for  ^  •„+!  s-nd 

Q.   -,   as  before  from  (3.5-7);   (ii)  we  correct   C^   ^   using 

(3.8-9)  which  requires  the  discharge  at  the  advanced  timestep. 

We  correct  C.   .   using  (3.8)  until  we  reach  the  right  boundary 

or  until  C,  ^-,  =  C.    to  within  the  accuracy  of  the  calculation, 
jn+1     jn  -^ 

The  correction  to   C .  ^,   was  carried  out  far  to  the  right  because 

jn+1  ^ 

numerical  experiments  showed  that  an  instability  appeared  at 

that  value  of  x,   which  was  the  last  corrected   C. 

J  jn+1 

obtained  by  using  (3.8).   When  the  correction  was  carried  far 
to  the  right,  no  instability  appeared, 
U.   Some  Numerical  Results. 

A  series  of  numerical  calculations  were  performed  on  an 

-  13  - 


1/2 


IBM  70k-   using  the  difference  scheme  described  above.   In  these 

calculations  fluid  was  Introduced  on  the  left  hand  boundary 

while  the  initial  flow  of  fluid  was  constant.   The  discharge 

Q.^   on  the  left  hand  boundary  was  of  the  form 
On  "^ 

(i;.10)  Q.Q^   =  bf-  gcos  inu^n^t)  "^  <  d^ 

the  initial  conditions  were 

(i;.2)  Q,.Q   =  b:-  S 

and 

(k.3)  c^-o  =  1 

and  the  right  boundary  was  taken  at  infinity. 

It  was  the  intention  of  these  calculations  to  see  if  fluid 
introduced  on  the  left  boundary  with  a  sufficiently  large  dis- 
charge  Q^   would  cause  a  wave  to  travel  down  the  channel  and 
^     On 

break.   In  all  our  calculations  we  did  not  find  any  evidence  of 
breaking!  even  when  the  free  stream  discharge  was  negative  (l,e, 
fluid  traveling  toward  the  left)  we  did  not  obtain  breaking.   In 
each  of  the  cases  breaking  would  have  occurred  long  before  our 
calculations  were  completed  according  to  the  usual  shallow 
water  theory  (cof.  Stoker  [l],  pp.  352-357)» 

In  Figure  2  we  see  the  build  up  of  C  .    and  its  development 
in  time  for  the  case  in  which  ^+ S  =  1.5   and  bC-)^=   1  , 
As  the  wave  of  fluid  traveled  down  stream  it  spread.   After  the 
wave  developed,  the  half  way  point  between  the  maximiim  height 
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and  the  free  stream  height  traveled  at  a  constant  speed   CT^  . 
The  computed  speed  CT^      agreed  very  closely  with  the  speed  of 
propagation  o""  given  by  titoker  [of.  Ref.  1,  pp.  321-323 J  for  a 
bore  with  height  equal  to  the  maximum  height  in  back  and  free 
stream  height  in  front.   In  Figure  3  we  have  plotted  the  wave 
profiles  against   x  -  d^t  ,   These  profiles  show  a  decided 
spread  very  similar  to  that  found  by  Gardner  and  Morikawa  L3] 
for  their  equation;  i.e.   c(x,  t)  =  fC'^^y.^  )  .   No  solutions  of 
this  type  have  been  found  analytically  for  equations  (2.5-6), 

In  Figure  U.   we  show  the  profiles  for  a  wave  arising  from 
a  large  discharge   Q^   at  the  left.   Again  no  breaking  was 
observed. 

In  Figure  5  we  show  the  profiles  for  a  wave  arising  from  a 
large  Inflow  discharge   Q„n  at  the  left  when  the  free  stream 
condition  corresponds  to  fluid  flowing  to  the  left.   We 
observed  no  breaking  and  the  wave  continued  to  grow  in  height. 
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